#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <list>
#include <cmath>

#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef double prec_t;

typedef TunaArray<prec_t, 2>::tiny floatTinyArray_t;
typedef TunaArray<prec_t, 2>::huge ScalarField2D;
typedef TunaArray<prec_t, 1>::huge ScalarField1D;

void writeToFile_DX(list<floatTinyArray_t> &, int, string);
void writeToFile_GNU(list<floatTinyArray_t> &, int, string);
void writeToFile_VTK(const ScalarField2D&, std::string, std::string);
void readFromFile_DX_2D(ScalarField2D&, ScalarField2D&, 
			int, string, int w = 5); 
void writeToFile_GNU(list<prec_t> &, int, string);

timer crono;

int main () 
{
    prec_t length_x, length_y, dt, left_extreme, right_extreme, y0, tol, dist;
    prec_t s0, cD;
    int num_nodes_x, num_nodes_y, Npoints, frec, Nstart, Nend, jump, Nxf, Nyf;

    prec_t etime = 0;
// -----> Reading data from the input file
    std::ifstream input_cin ("input01");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> y0
	      >> left_extreme
	      >> right_extreme
	      >> Npoints
	      >> Nstart
	      >> Nend
	      >> jump
	      >> dt
	      >> s0
	      >> cD
	      >> tol
	      >> Nxf
	      >> Nyf
	      >> frec;
    input_cin.close();

    // THE TIME STEP MUST BE CALCULATED ACCORDING TO THE JUMP
    dt = dt * jump;

    /*
    cout << "\n +-----+ \n" 
	 << length_x << "\t" << length_y << "\t" << num_nodes_x << "\t" << num_nodes_y << "\n"
	 << left_extreme << "\t" << right_extreme << "\t" << Npoints << "\n"
	 << Nstart << "\t" << Nend << "\t"<< dt << "\n"
	 << tol << "\t" << frec 
	 << "\n +-----+ \n";
	 /**/

// -----> Mesh to read the velocity
    StructuredMesh< Uniform<prec_t, 2> > mesh(length_x, num_nodes_x, 
					      length_y, num_nodes_y);
    floatTinyArray_t deltas = mesh.getDeltas();
    //    mesh.print();
    
// -----> List for storing the line points. 
    list<floatTinyArray_t> linec, linec0;
    list<floatTinyArray_t>::iterator iter, iter0;
    floatTinyArray_t aux; 
    floatTinyArray_t x_i_n;     // Coordinates of x(i) at step n
    floatTinyArray_t x_i_np1;   // Coordinates of x(i) at step n+1
    floatTinyArray_t x_im1_n;   // Coordinates of x(i-1) at step n
    floatTinyArray_t x_im1_np1; // Coordinates of x(i-1) at step n+1
    floatTinyArray_t x_i_0, x_im1_0;

//   Strip geometry
//  ---------------------------------------
//        deltaL0
//        |-----|
//        v     v
//  +-----------------------+   <-+
//  +-----+-----+-----+-----+     | s0
//  +-----------------------+   <-+
//
//  ---------------------------------------
    prec_t deltaL0 = (right_extreme - left_extreme) / (Npoints - 1);
    prec_t deltaL0R = deltaL0;
    prec_t line_length = 0.0;

    //    ofstream flen("/data/tmpu/DataMix2D_01L/lineLength"); 
    //    ofstream fPDF("/data/tmpu/DataMix2D_01L/PDF"); 

    cout << "\n +----- Line construction -----+ ... ";

    list<prec_t> s, tao, rho, delta;
    list<prec_t>::iterator iS, iT, iR, iD;

    for(int i = 0; i < Npoints; i++) {
	aux(0) = left_extreme + i * deltaL0;
	aux(1) = y0;
	linec.push_back(aux);
	linec0.push_back(aux);

	// We need S and Tao for each point (interval)
	s.push_back(s0);
	tao.push_back(0.0);
	rho.push_back(1.0);
	delta.push_back(deltaL0);
    }

    cout << "\n Sizes: line = " << linec.size() 
	 << "; s = " << s.size() 
	 << "; tao = " << tao.size()
	 << "; rho = " << rho.size() << "\n" << flush;

    cout << " +----- DONE -----+\n";

    ////////////////////// 
    /* *
    iter = linec.begin(); iS = s.begin(); iT = tao.begin();
    iR = rho.begin(); iD = delta.begin();
    for(int i = 0; i < Npoints; i++) {
      cout << i << " " << (*iter)(0) << " " << (*iter)(1) << " "
	   << *iS << " " << *iT << " " << *iR << " " << *iD << endl;
      iter++; iS++; iT++; iR++; iD++;
    }
    int pausa; cin >> pausa;
    /* */
    ////////////////////// 

    writeToFile_DX(linec, 0, "/data/tmpu/DataMix2D_01L/line.");

    cout << "\n +---------------------------------------+"
	 << "\n +      Beginning the line tracking      |"
	 << "\n +---------------------------------------+\n\n";

// -----> Velocity field at time t
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());

// -----> Velocity field at time t + dt
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());

//  Particle tracking: 1) Linear interpolation of velocity, 
//                     2) Runge-Kutta fourth order integrator
//                     3) 2D, prec_t = precision.   
    ParticleTracking< RungeKutta4< Linear< prec_t, 2 > > > particle;

    ScalarField2D c(Nxf, Nyf);
    prec_t sigmax, sigmay, normSig, coordx, coordy, midx, midy;
    prec_t aux1, aux2, auxDelta1, auxDelta2;
    prec_t dxf = length_x / (Nxf - 1);
    prec_t dyf = length_y / (Nyf - 1);
    prec_t lambda;

    int Nnew, Ntemp;
    prec_t d0, d1, frac;

    cout << " +----- Time series -----+\n"
	 << " | Nstar = " << Nstart << "\t Nend = " << Nend 
	 << "\t jump = " << jump 
	 << "\n +-----------------------+ ";

// +------------------------------------------+
// -----> Main loop
// +------------------------------------------+
    int contador = 0;
    double timeTrack, totalTimeTrack = 0;
    double timeDiff, totalTimeDiff = 0;
    double timeFileRead, totalTimeFileRead = 0;
    double desvTrack[(Nend - Nstart) / jump];
    double desvDiff[(Nend - Nstart) / jump];
    double desvFile[(Nend - Nstart) / jump];
    for(int step = Nstart; step < Nend; step += jump, contador++) {
      
      cout << "\n\n ----> Step = " << step ;
      /* *
      cout << "\t Npoints = " << linec.size() << flush;
      cout << "\n Sizes: line = " << linec.size() 
	   << "; s = " << s.size() 
	   << "; tao = " << tao.size()
	   << "; rho = " << rho.size() << "\n" << flush;
      /* */


      cout << "\n   +----- Vec file read ...  "; 
      crono.tic(); // Read file
// -----> Read the velocity at time t from a file
      readFromFile_DX_2D(u, v, step, "/data/tmpu/Data4/velc.");
// -----> Read the velocity at time t + dt from a file
      readFromFile_DX_2D(un, vn, step + jump, "/data/tmpu/Data4/velc.");
      timeFileRead = crono.toc();
      totalTimeFileRead += timeFileRead;
      cout << " Elapsed time = " << timeFileRead << flush ;

      cout << "\n   +----- Tracking ... " << flush;
      crono.tic(); // Tracking
// -----> First point on the line
      iter = linec.begin(); // Pointer to initial point in the line x(i-1)^n
      x_im1_n = *iter;      // Contents of initial point in the line x(i-1)^n

// -----> Track first point to the next position x(i-1)^(n+1) 
      *iter = particle.track(*iter, deltas, u, v, un, vn, dt);
      x_im1_np1 = *iter;
      iter++;  // Advance iter to the second point on the line

      line_length = 0.0;
      deltaL0 = 0.0;

      iS = s.begin(); iT = tao.begin(); iD = delta.begin();
      iS++; iT++; iD++; 

// -----> Track the rest of particles on the line
      for( ; 
	   iter != linec.end(); 
	   ++iter, ++iS, ++iT, ++iD)  {

	x_i_n = *iter; // Current point  x(i)^n
	x_i_np1 = particle.track(*iter, deltas, u, v, un, vn, dt); 

/* */
// -----> Calculate the distance between consecutive tracked points 
//        x(i-1)^(n+1) and x(i)^(n+1)
	dist = sqrt( ( x_i_np1(0) - x_im1_np1(0) ) * 
		     ( x_i_np1(0) - x_im1_np1(0) ) +
		     ( x_i_np1(1) - x_im1_np1(1) ) * 
		     ( x_i_np1(1) - x_im1_np1(1) ) );	

// -----> Add points if dist exceeds the tolerance, 
//        adds a point between points: x(i)^n and x(i-1)^n
/* */
	if ( dist > tol ) {

	  aux(0) = ( x_i_n(0) + x_im1_n(0) ) * 0.5;
	  aux(1) = ( x_i_n(1) + x_im1_n(1) ) * 0.5;

// -----> Insert one point, before the actual position of iter
	  linec.insert(iter,aux); 
// -----> Rewind the iterator 2 positions, (the for loop will advance one pos.)
	  iter--; iter--; 

// -----> Striation thickness will keep its value for the new two intervals.
	  s.insert(iS, *iS); 
// -----> Rewind iS two positions, to the previously inserted striation thickness
	  iS--;iS--; 

// -----> Dimensionless time will keep its value for the new two intervals.
	  tao.insert(iT, *iT);
// -----> Rewind iT two positions, to the previously inserted tao. 
	  iT--;iT--; 

// -----> Insert the half delta for the new inserted point
	  delta.insert(iD, dist * 0.5);
// -----> Rewind iD two positions, to the previously inserted delta. 
	  iD--; iD--;

	} else { // If dist does not exceeds the tolerance, 		
/* */
// -----> Re-state the value of x(i-1)^n
	  x_im1_n   = x_i_n;   
// -----> Re-state the value of x(i-1)^(n+1)
	  x_im1_np1 = x_i_np1; 
// -----> Modify the new position of x(i) in the list ( x(i)^(n+1) )
	  *iter     = x_i_np1; 	  
// -----> Striation thickness: s_i = s0 * delta_0 / delta_i
	  *iS =  s0 * deltaL0R / dist; 
// -----> Calculate the dimensionless time
	  *iT = *iT + dt * cD / ( (*iS) * (*iS) ); 

	  deltaL0 += *iS * sqrt(1 + 4 * (*iT)); 
// Accumulate the length of the strip	 	  		  
	  line_length += dist; 
/* */
        }
/* */
      
      }

      deltaL0 /= linec.size();
      /* *
      cout << "\n\t DeltaL0 = " << deltaL0 
	   << "\n\t Length = " << line_length  << "\n\t ... "<< flush; 
      /* */
      //      flen << step << "\t" << line_length << endl;
      timeTrack = crono.toc();
      cout << " Elapsed time = " << timeTrack << flush ;
      totalTimeTrack += timeTrack;
      
      /* *
// ------> Reinterpolation


      deltaL0 = deltaL0R * 0.3;


//      cout << "\n\t DeltaL0R = " << deltaL0R << "\t DeltaL0 = " << deltaL0 
//	   << "\t length = " << line_length << flush; 


      if ( !(step % 500) ) {
	cout << "\n   +----- Reinterpolation ... " << flush;
	crono.tic();
	
	Nnew = int(line_length / deltaL0);
	frac = (line_length - Nnew * deltaL0) / deltaL0;      
	if ( frac > 0.5 ) Nnew++;
	deltaL0 = line_length / Nnew;
	
	//	cout << "\n\t Nnew = " << Nnew << "\t delta =" << deltaL0 << flush; 
	
	iter = linec.begin(); // Pointer to initial point in the line x(i-1)^n
	x_im1_n = *iter;      // Contents of initial point in the line x(i-1)^n
	iter++;  // Advance iter to the second point in the line
	
	iS = s.begin(); iT = tao.begin(); iD = delta.begin();
	iS++; iT++; iD++; 
	
	d0 = 0.0; d1 = 0.0; auxDelta2 = 0.0;
	for( ; 
	     iter != linec.end(); 
	     ++iter, ++iS, ++iT, ++iD)  
	  {
	  
	    x_i_n = *iter;
	    
	    auxDelta1 = *iD - d1;
	    Ntemp = int( auxDelta1 / deltaL0 );
	    
	    cout << "\n Ntemp = " << Ntemp << "\t auxDelta = " << auxDelta1 << "\t deltaL0 = " << deltaL0 << flush;
	    cout << "\n d0 = " << d0 << "\t d1 = " << d1 << "\t x(i-1) = " << x_im1_n(0) << flush;
	    
	    if (d1 > 0) {
	      
	        cout << "\n d1 > 0"; 
		lambda = d1 / (*iD);
		x_i_np1(0) = lambda * x_i_n(0) + (1 - lambda) * x_im1_n(0);
		x_i_np1(1) = lambda * x_i_n(1) + (1 - lambda) * x_im1_n(1);
		iter--; iS--; iT--; iD--;
		*iter = x_i_np1;
		*iS = 0;
		*iT = 0;
		*iD = deltaL0;
		iter++; iS++; iT++; iD++;
				
		cout << "\n x(i) = " << x_i_n(0) << "\t x(i-1) = " << x_im1_n(0) << flush;		
		cout << "\n Int x(i) = " << x_i_np1(0) << "\n" << flush;  
		
		for(int j = 1; j <= Ntemp; ++j) {
		  lambda = (d1 + deltaL0 * j) / (*iD);
		  x_i_np1(0) = lambda * x_i_n(0) + (1 - lambda) * x_im1_n(0);
		  x_i_np1(1) = lambda * x_i_n(1) + (1 - lambda) * x_im1_n(1);	
		  
		  //		cout << "\n " << j << " x(i) = " << x_i_np1(0) << flush;  
		  
		  linec.insert(iter, x_i_np1);
		  delta.insert(iD, deltaL0);
		  s.insert(iS, 0);
		  tao.insert(iT, 0);
		} 

	    } else {
	      
	        cout << "\n d1 <== 0 "; 
	      
		lambda = deltaL0 / (*iD);	       
		
		for(int j = 1; j <= Ntemp; ++j) {
		  x_i_np1(0) = lambda*j * x_i_n(0) + (1 - lambda*j) * x_im1_n(0);
		  x_i_np1(1) = lambda*j * x_i_n(1) + (1 - lambda*j) * x_im1_n(1);	  
		  
		  cout << "\n x(i) = " << x_i_n(0) << "\t x(i-1) = " << x_im1_n(0) << flush;		
		  cout << "\n Int x(i) = " << x_i_np1(0) << "\n" << flush;  
		  
		  linec.insert(iter, x_i_np1);
		  delta.insert(iD, deltaL0);
		  s.insert(iS, 0);
		tao.insert(iT, 0);
		}
	    }
	    d0 = auxDelta1 - deltaL0 * Ntemp;
	    d1 = deltaL0 - d0;	  
	    x_im1_n = *iter;
	  
	  
	  }
      }
      //cout << endl;
      cout << " Elapsed time = " << crono.toc() << flush ;

      /* */
      if( !(step % frec) ) {
	cout << "\n   +----- Writing line." << step <<" file ... " << flush ;
	crono.tic();
	writeToFile_DX(linec, step, "/data/tmpu/DataMix2D_01L/line.");
	cout << " Elapsed time = " << crono.toc() << flush ;
      }
      /* */
      cout << "\n   +----- Diffusion calculation ... " << flush;


      /* *
//---> This is for the PDF <----------------------------------------------------------------
      iter0 = linec0.begin(); 
      *iter0 = particle.track(*iter0, deltas, u, v, un, vn, dt);
      x_im1_0 = *iter0;
      iter0++;  // Advance iter to the second point in the line
      for( iR = rho.begin(); 
	   iter0 != linec0.end(); 
	   ++iR, ++iter0)  {
	x_i_0 = particle.track(*iter0, deltas, u, v, un, vn, dt); 	  
	*iR = sqrt( ( x_i_0(0) - x_im1_0(0) ) * 
		      ( x_i_0(0) - x_im1_0(0) ) +
		      ( x_i_0(1) - x_im1_0(1) ) * 
		      ( x_i_0(1) - x_im1_0(1) ) ) / deltaL0R;
	x_im1_0 = x_i_0; 
	*iter0  = x_i_0; 	
      }

      if ( !(step % 100) ) {
	rho.sort();
	cout << "\n Rho = " << flush;
	for( iR = rho.begin(); iR != rho.end(); ++iR) 
	  cout << *iR << "\t" << flush;
	cout << "\n" << flush;
	writeToFile_GNU(rho, step, "/data/tmpu/DataMix2D_01L/PDF.");
      }
//------------------------------------------------------------------------------------------ 

/* *

// ------> Calculate the Delta L for the reinterpolation
      crono.tic();
      deltaL0 = 0.0;

      iS = s.begin(), iT = tao.begin();
      iS++; iT++;

      for(; iS != s.end(); ++iS, ++iT) 
	deltaL0 += *iS * sqrt(1 + 4 * (*iT)); 
      deltaL0 /= linec.size();

      cout << " DeltaL = " << deltaL0 << " " << flush; 

      
      /* */

// ------> Reconstruct the concentration field on a fine grid 

      crono.tic(); // Diffusion
      c = 0;

      /*
      iT = tao.begin(); iT++;
      float AUXTAO = *iT;
      iT = tao.begin(); *iT = AUXTAO;
      */

      for(int i = 0; i < Nxf; ++i)
	for(int j = 0; j < Nyf; ++j) {
	  iter = linec.begin();  
	  x_im1_n = *iter; 
	  iter++;
	  iS = s.begin(), iT = tao.begin();
	  iS++; iT++;
	  for( ; iter != linec.end(); ++iter, ++iS, ++iT) {
	    x_i_n = *iter;
	    sigmax = x_i_n(0) - x_im1_n(0); 
	    sigmay = x_i_n(1) - x_im1_n(1);
	    normSig = sqrt ( sigmax * sigmax + sigmay * sigmay);
	    sigmax /= normSig;
	    sigmay /= normSig;
	    coordx = dxf * i;
	    coordy = dyf * j;
	    midx = 0.5 * ( x_i_n(0) + x_im1_n(0) ); 
	    midy = 0.5 * ( x_i_n(1) + x_im1_n(1) );	  
	    aux1 = (coordx - midx) * sigmax + (coordy - midy) * sigmay;
	    aux2 = (coordx - midx) * (-sigmay) + (coordy - midy) * sigmax;
	    aux1 = aux1 * aux1 / (deltaL0 * deltaL0);
	    aux2 = aux2 * aux2 / ( (*iS) * (*iS) * (1.0 + 4.0 * (*iT)) );
	    c(i,j) += exp(-aux1-aux2) / (1.7726 * sqrt( 1.0 + 4.0 * (*iT)));
	    x_im1_n = *iter;
	  }
	}

      timeDiff = crono.toc();
      cout << " Elapsed time = " << timeDiff << flush ;
      totalTimeDiff += timeDiff;
      
      auxDelta1 = max(c);
      c = c / auxDelta1;

      /* */
      cout << "\n   +----- Writing conc." << step << " file ..." << flush;
      crono.tic();
      InOut::writeToFile_DX(c, step, "/data/tmpu/DataMix2D_01L/conc.", dxf, dyf);  
      cout << " Elapsed time = " << crono.toc() << flush;

      cout << "\n max = " << auxDelta1 << endl;

      desvTrack[contador] = timeTrack;
      desvDiff[contador] = timeDiff;
      desvFile[contador] = timeFileRead;
      
      /* */

    }

// +------------------------------------------+
// -----> End of main loop
// +------------------------------------------+

//    flen.close();
//    fPDF.close();


    double AVTrack = totalTimeTrack / contador;
    double AVDiff  = totalTimeDiff / contador;
    double AVFile  = totalTimeFileRead / contador;

    double DSTrack = 0;
    double DSDiff = 0;
    double DSFile = 0;
    for(int i = 0; i < contador; i++) {
      DSTrack += (AVTrack - desvTrack[i])*(AVTrack - desvTrack[i]);
      DSDiff += (AVDiff - desvDiff[i])*(AVDiff - desvDiff[i]);
      DSFile += (AVFile - desvFile[i])*(AVFile - desvFile[i]);
    }
    DSTrack = sqrt(DSTrack / (contador - 1) );
    DSDiff = sqrt(DSDiff / (contador - 1) );
    DSFile = sqrt(DSFile / (contador - 1) );

    cout << "\n +------------------------------------------+ "
	 << "\n | Total time track = " << totalTimeTrack
	 << "\n | Total time diff = " << totalTimeDiff
	 << "\n | Total time file read = " << totalTimeFileRead
	 << "\n | Average track = " << AVTrack
	 << "\n | Average diff = " << AVDiff
	 << "\n | Average read = " << AVFile
	 << "\n | Desviacion track = " << DSTrack
	 << "\n | Desviacion diff = " << DSDiff
	 << "\n | Desviacion file = " << DSFile
	 << "\n +------------------------------------------+ \n";

    cout << "\n End ... \n\n";

    return 0;
}

void writeToFile_DX(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    //    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 2 items "
	 << Npoints << " data follows\n";
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "   " << (*iter)(1) << "   ";
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << Npoints - 1 << " data follows\n";
    for(int i = 0; i < Npoints - 1; ++i)
	file << i << "   " << i + 1 << "   ";
    file << endl;

    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";
    file << "object 3 class array type float rank 0 items "
	 << Npoints << " data follows\n";
    for(int i = 0; i < Npoints; ++i)
	file << 1 << "\t";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();

}

void writeToFile_GNU(list<prec_t> &lista, int n, string filename)
{
    int Npoints = lista.size();
    list<prec_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    //    cout << "writing " << filename << endl;

    int N = 10;
    prec_t arreglo[N];
    for(int i = 0; i < N; ++i) arreglo[i] = 0.0;

    iter = lista.begin();
    prec_t min = *iter;
    iter = lista.end(); --iter;
    prec_t max = *iter;
    prec_t delta_rho = (max - min) / N;

    for( iter = lista.begin(); iter != lista.end(); ++iter)
      for(int i = 0; i < N; ++i)
	if ( *iter >= (min + i * delta_rho) &&
	     *iter < min + (i+1) * delta_rho)
	  arreglo[i] += 1;
    
    for(int i = 0; i < N; ++i)
      file << (min + i * delta_rho) * 0.5 << "\t" << arreglo[i] << "\n";
    file << endl;

    /*
    int i = 1;
    for(iter = lista.begin(); iter != lista.end(); ++iter, ++i)
      file << i << "\t" << *iter << "\n";
    file << endl;
    */

    file.close();

}

void writeToFile_GNU(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "\t" << (*iter)(1) << "\n";
    file << endl;

    file.close();

}

void writeToFile_VTK(const ScalarField2D& data, 
		     std::string name, std::string variable)
{
    int Nvtk_x = data.rows();
    int Nvtk_y = data.cols();
    int Nvtk_z = 1;
    prec_t dx = 1. / (Nvtk_x - 1);
    prec_t dy = 1. / (Nvtk_y - 1);
    prec_t dz = 1.;
    int Ntotal = Nvtk_x * Nvtk_y * Nvtk_z;

    std::ofstream file (name.c_str());        
    file << "# vtk DataFile Version 2.0 \n"
	 << "Lid_driven streamfuntion-vorticity \n"
	 << "ASCII \n"
	 << "DATASET RECTILINEAR_GRID \n" 
	 << "DIMENSIONS " 
	 << Nvtk_x << " " << Nvtk_y << " " << Nvtk_z;

    file << "\nX_COORDINATES " << Nvtk_x << " float \n";
    for(int i = 0; i < Nvtk_x; ++i) {
	file << i * dx << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nY_COORDINATES " << Nvtk_y << " float \n"; 
    for(int i = 0; i < Nvtk_y; ++i) {
	file << i * dy << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nZ_COORDINATES " << Nvtk_z << " float \n"; 
    for(int i = 0; i < Nvtk_z; ++i) {
	file << i * dz << " ";
	if ( !(i % 6) ) file << "\n";
    }
	
    file << "\nPOINT_DATA " << Ntotal
	 << "\nSCALARS " << variable << " float"
	 << "\nLOOKUP_TABLE default \n";
    for(int k = 1; k <= Nvtk_z ; ++k) {
	for(int j = 1; j <= Nvtk_y ; ++j) {
	    for(int i = 1; i <= Nvtk_x; ++i) {
		file << data(i,j) << " ";
		if ( !(i % 6) ) file << "\n";
	    }
	}
    }
        
    file.close();

}

void readFromFile_DX_2D(ScalarField2D& scf1, ScalarField2D& scf2, 
			int number, string filename, int w) 
{   
  //  int w = 5;
  ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;
    ifstream file (filename.c_str());

    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    } else {
      //      cout << "\n Opening file \" " << filename << "\"" << endl;
    }
    
    char line[256];
/// 
/// by the time, just jump the first 8 lines
///	
    for(int ii = 0; ii < 8; ++ii) {
	file.getline(line,256);
	//	cout << " ii = " << ii << "\t" << line << endl;
    }

    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
///
/// Here, scf1 and scf2 are supposed to be of type ScalarField2D
///
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j)
	  file >> scf1(i,j) >> scf2(i,j) ;
    file.close();
}


